#pragma OPENCL EXTENSION cl_khr_fp64: enable

typedef struct {
	FLOAT real;
	FLOAT imag;
} Complex;

//#define Complex double2
//#define real x
//#define imag y

inline Complex complex_number(FLOAT real, FLOAT imag) {
	Complex c;
	c.real = real;
	c.imag = imag;
	return c;
}

inline FLOAT real(Complex c) {
	return c.real;
}

inline FLOAT imag(Complex c) {
	return c.imag;
}

inline Complex complex_add(Complex a, Complex b) {
	Complex c;
	c.real = a.real + b.real;
	c.imag = a.imag + b.imag;
	return c;
}

inline Complex complex_sub(Complex a, Complex b) {
	Complex c;
	c.real = a.real - b.real;
	c.imag = a.imag - b.imag;
	return c;
}

inline Complex complex_mul(Complex a, Complex b) {
	Complex c;
	c.real = a.real * b.real - a.imag * b.imag;
	c.imag = a.real * b.imag + a.imag * b.real;
	//c.real = b.real * b.real - a.imag * a.imag;
	//c.imag = a.real * b.imag + a.imag * b.real;
	return c;
}

inline Complex complex_div(Complex a, Complex b) {
	Complex c;
	c.real = (a.real * b.real + a.imag * b.imag) / (b.real * b.real + b.imag * b.imag);
	c.imag = (a.imag * b.real - a.real * b.imag) / (b.real * b.real + b.imag * b.imag);
	return c;
}

// FIXME: this can only raise complex numbers up to an integer power.
// in future we may want to raise to a non-integer real number and/or a complex number
/*
inline Complex complex_pow(Complex a, int b) {
	Complex c = a;
	for (int i = 1; i < b; i++) {
		c = complex_mul(c, a);
	}
	return c;
}
*/
inline Complex complex_pow(Complex a, Complex b) {
	Complex c = a;
	for (int i = 1; i < b.real; i++) {
		c = complex_mul(c, a);
	}
	return c;
}

inline FLOAT complex_mod(Complex a){
    return (sqrt(a.real * a.real + a.imag * a.imag));
}

inline FLOAT complex_arg(Complex a){
	if (a.real > 0) {
		return atan(a.imag / a.real);
	} else if (a.real < 0 && a.imag >= 0) {
		return atan(a.imag / a.real) + M_PI;
	} else if (a.real < 0 && a.imag < 0) {
		return atan(a.imag / a.real) - M_PI;
	} else if (a.real == 0 && a.imag > 0) {
		return M_PI / 2;
	} else if (a.real == 0 && a.imag < 0) {
		return -M_PI / 2;
	} else {
		return 0;
	}
}

inline Complex complex_sqrt(Complex a) {
	Complex c;
	//(cfloat)( sqrt(cmod(a)) * cos(carg(a)/2),  sqrt(cmod(a)) * sin(carg(a)/2));
	c.real = sqrt(complex_mod(a)) * cos(complex_arg(a) / 2);
	c.imag = sqrt(complex_mod(a)) * sin(complex_arg(a) / 2);
	return c;
}

inline Complex complex_conj(Complex c) {
	c.imag *= -1.;
	return c;
}

__kernel void fractal(int w, FLOAT scale, FLOAT topLeftX, FLOAT topLeftY, int limit, __global FLOAT *buffer) {
	int i;
	//const FLOAT escapeRadius = 2;
	int x = get_global_id(0);
	int y = get_global_id(1);
	Complex z;
	Complex c;
	z.real = 0; z.imag = 0;
	c.real = (x * scale) + topLeftX;
	c.imag = (y * scale) + topLeftY;
	for (i = 0; (i < limit) && ((z.real * z.real + z.imag * z.imag) < 4); i++) {
		// LOOP BODY
	}
	buffer[y * w + x] = ((FLOAT)i) / limit;
}

// NB: "FLOAT" is just a macro, usually defined to "double", but could be "float" if desired.
__kernel void julia(int w, FLOAT scale, FLOAT topLeftX, FLOAT topLeftY, FLOAT cx, FLOAT cy, int limit, __global FLOAT *buffer) {
	int i;
	int x = get_global_id(0);
	int y = get_global_id(1);
	Complex z;
	Complex c;
	z.real = (x * scale) + topLeftX;
	z.imag = (y * scale) + topLeftY;
	c.real = cx; c.imag = cy;
	
	for (i = 0; (i < limit) && ((z.real * z.real + z.imag * z.imag) < 4); i++) {
		// LOOP BODY
	}
	buffer[y * w + x] = ((FLOAT)i) / limit;
}
